/* 
* 
* 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 "esriheader.h"
#include "Ridge.h"
#include "Erreur.h"
#include <set>
#include <map>
#include <algorithm>

using namespace std;


// Definition : a pixel belongs to ridge set R(a,k) iff { it belongs to the boundary of a catchment draining at least "a"
//		and it is located at a distance of at least k*sqrt(a) from the outlet of this catchment }

int GetUpstream2 ( FlowDirGrid * F , int * _cell , int * _nup , int ** _UPSTREAM , set<int> & Neighbors )
{
    int Dr,Dc,nr,nc,Drmin,Drmax,Dcmin,Dcmax,irec,nup;
    
    *_UPSTREAM = NULL;
    nup = 0;
    
    (void) rec2RC ( &(F->Header) , _cell , & nr , & nc );
    (void) GetNeighborhood ( &(F->Header) , & nr , & nc , & Drmin , & Drmax , & Dcmin , & Dcmax );

    Neighbors.clear();

    for(Dr=Drmin;Dr<=Drmax;Dr++)
    {
        for(Dc=Dcmin;Dc<=Dcmax;Dc++)
        {
			if ( Dr !=0 | Dc !=0 )
			{
				irec = (nr+Dr)*(F->Header.ncols)+(nc+Dc);
				Neighbors.insert(irec);

				if ( F->Data[nr+Dr][nc+Dc] == F->Compass[1-Dr][1-Dc] )
				{
					if ( !( (* _UPSTREAM) = (int *) realloc( (* _UPSTREAM) ,(nup+1)*sizeof(int) ) ) )
					{
						SetErreur("GetUpstream : could not reallocate integer pointer\n");
						return(-1);
					}

					(*_UPSTREAM)[nup] = irec;
					nup += 1;
				}	
            }
        }
    }

    *_nup = nup;

    return(0);

}

int GetNeighborhood2 ( esriheader * Header , int * _cell , set<int> & Neighbors )
{
    int Dr,Dc,nr,nc,Drmin,Drmax,Dcmin,Dcmax,irec,nup;

    (void) rec2RC ( Header , _cell , & nr , & nc );

    if( nr == 0 )
        Drmin = 0;
    else
        Drmin = -1;

    if( nr == (Header->nrows) -1 )
        Drmax = 0;
    else
        Drmax = 1;

    if( nc == 0 )
        Dcmin = 0;
    else
        Dcmin = -1;

    if( nc == (Header->ncols) -1 )
        Dcmax = 0;
    else
        Dcmax=1;

    Neighbors.clear();

    for(Dr=Drmin;Dr<=Drmax;Dr++)
    {
        for(Dc=Dcmin;Dc<=Dcmax;Dc++)
        {
			if ( Dr !=0 | Dc !=0 )
			{
				irec = (nr+Dr)*(Header->ncols)+(nc+Dc);
				Neighbors.insert(irec);
			}
        }
    }

    return(0);

}


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

int TestCast ( int * MYINTS , int n, float ** _MYFLOATS ) 
{
    int myint, i ;  
    float myfloat ;
    float * ptrfloat ;

    *_MYFLOATS = (float *) malloc ( n * sizeof(float) ) ;

    for(i=0;i<n;i++)
    {
            myint = MYINTS[i] ;
            ptrfloat = (float *) ( & myint ) ;
	    (*_MYFLOATS)[i] = *ptrfloat ;   
    }
   
    return(0);
}

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


int TestEnsembles ( 	int * INNERBOUNDARY , int ni ,
			int * OUTERBOUNDARY , int no ,
			int * INNERSUBBOUNDARY , int nis ,
			int * OUTERSUBBOUNDARY , int nos ,
			int ncols , int nrows ,
			int ** _INNERBOUNDARY , int * _ni ,
			int ** _OUTERBOUNDARY , int * _no )
{

	set<int> InnerBoundary, OuterBoundary, InnerSubBoundary, OuterSubBoundary ;
	set<int> TmpSet1, TmpSet2 ;
	set<int>::iterator isi, isi2 ;
	int i,nr,nc,irec,irec2,Drmin,Drmax,Dcmin,Dcmax,Dr,Dc ;

	InnerBoundary.clear();

	for(i=0;i<ni;i++)
		InnerBoundary.insert(INNERBOUNDARY[i]);

	OuterBoundary.clear();

	for(i=0;i<no;i++)
		OuterBoundary.insert(OUTERBOUNDARY[i]);
		 
	InnerSubBoundary.clear();

	for(i=0;i<nis;i++)
		InnerSubBoundary.insert(INNERSUBBOUNDARY[i]);

	OuterSubBoundary.clear();

	for(i=0;i<nos;i++)
		OuterSubBoundary.insert(OUTERSUBBOUNDARY[i]);


	// Test ensemble operation 

		TmpSet1.clear();
		std::set_intersection (OuterBoundary.begin(),OuterBoundary.end(), InnerSubBoundary.begin(),InnerSubBoundary.end(), 
						std::insert_iterator< std::set<int> >( TmpSet1, TmpSet1.begin() ) );
		std::set_intersection (OuterSubBoundary.begin(),OuterSubBoundary.end(), InnerBoundary.begin(),InnerBoundary.end(), 
						std::insert_iterator< std::set<int> >( TmpSet1, TmpSet1.end() ) );

		// Merge the two outer boundaries and remove TmpSet1 from the union

		OuterBoundary.insert(OuterSubBoundary.begin(),OuterSubBoundary.end());	
//		for(isi=OuterSubBoundary.begin();isi!=OuterSubBoundary.end();isi++)
//			OuterBoundary.insert((*isi));

		TmpSet2.clear();
		std::set_difference (OuterBoundary.begin(),OuterBoundary.end(), TmpSet1.begin(), TmpSet1.end(),
						std::insert_iterator< std::set<int> >( TmpSet2, TmpSet2.begin() ) );

		OuterBoundary.clear();
		OuterBoundary.insert(TmpSet2.begin(),TmpSet2.end());
//		for(isi=TmpSet2.begin();isi!=TmpSet2.end();isi++)
//			OuterBoundary.insert((*isi));

		// Do the same for inner boundaries

		InnerBoundary.insert(InnerSubBoundary.begin(),InnerSubBoundary.end());
//		for(isi=InnerSubBoundary.begin();isi!=InnerSubBoundary.end();isi++)
//			InnerBoundary.insert((*isi));

		TmpSet2.clear();
		std::set_difference (InnerBoundary.begin(),InnerBoundary.end(), TmpSet1.begin(), TmpSet1.end(),
						std::insert_iterator< std::set<int> >( TmpSet2, TmpSet2.begin() ) );

		InnerBoundary.clear();
		InnerBoundary.insert(TmpSet2.begin(),TmpSet2.end());
//		for(isi=TmpSet2.begin();isi!=TmpSet2.end();isi++)
//			InnerBoundary.insert((*isi));

		// Some points in TmpSet1 may still belong to the new inner boundaries
		// 	test if they have a neighbor in the final OuterBoundary

		for (isi=TmpSet1.begin();isi!=TmpSet1.end();isi++)
		{
			irec = (*isi);
		    	nc = irec%ncols ;
		    	nr = (irec - nc)/ncols ;

//			SetErreur("Analyse du voisinage de %d (row %d, col %d)\n",irec,nr,nc);

		    	if( nr == 0 )
				Drmin = 0;
		    	else
				Drmin = -1;

			if( nr == (nrows-1))
				Drmax = 0;
			else
				Drmax = 1;

			if( nc == 0 )
				Dcmin = 0;
		    	else
				Dcmin = -1;

		    	if( nc == (ncols-1))
				Dcmax = 0;
		    	else
				Dcmax=1;

			for(Dr=Drmin;Dr<=Drmax;Dr++)
			{
				for(Dc=Dcmin;Dc<=Dcmax;Dc++)
				{
					if ( Dr !=0 | Dc !=0 )
					{
						irec2 = (nr+Dr)*ncols+(nc+Dc);
						isi2 = OuterBoundary.find(irec2);
//						SetErreur("   Analyse du voisin %d\n",irec2);

						if (isi2!=OuterBoundary.end() )
						{
							InnerBoundary.insert(irec);
//							SetErreur("       Le voisin %d est dans OuterBoundary\n",irec2);
						}
					}
				}
			}
		}


	// Store results
	
	(*_ni) = InnerBoundary.size();
	(*_INNERBOUNDARY) = (int *) malloc ( (*_ni) * sizeof(int) );

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

	(*_no) = OuterBoundary.size();
	(*_OUTERBOUNDARY) = (int *) malloc ( (*_no) * sizeof(int) );

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

	InnerBoundary.clear();
	OuterBoundary.clear();
	InnerSubBoundary.clear();
	OuterSubBoundary.clear();
	TmpSet1.clear();
	TmpSet2.clear();

	return(0);

} 

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

int RidgeMagnitudeArea_ ( FlowDirGrid * F , FloatRaster * D , float * _k , FloatRaster * M ,
			    int * _root , set<int> & InnerBoundary , set<int> & OuterBoundary , int * _nup ,
			    set<int> & TmpSet1 , set<int> & TmpSet2 ,
			    set<int>::iterator & isi , set<int>::iterator &  isi2, set<int>::iterator &  isi3 ,
			    int * _irec , int * _nr , int * _nc , float * _d0 , float * _d , float * _a )			
{
	set<int> InnerSubBoundary , OuterSubBoundary ; //, TmpSet1 , TmpSet2 ;
	int * UPSTREAM ;
//	float d0, d , a ;
	int i,nup1,nupi ;//irec,nr,nc;
//	set<int>::iterator isi, isi2, isi3;

	OuterBoundary.clear();
	InnerBoundary.clear();

	InnerBoundary.insert(*_root);
	(*_nup) = 1;

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

	// Extraction of upstream subcatchments
	
	for(i=0;i<nup1;i++)
    	{
		if ( RidgeMagnitudeArea_ ( F , D , _k , M , UPSTREAM+i , InnerSubBoundary , OuterSubBoundary , & nupi ,
					TmpSet1 , TmpSet2 , isi , isi2 , isi3 , _irec , _nr , _nc , _d0 , _d , _a ) )
		{
			SetErreur("RidgeMagnitude_ : recursive call failed\n");
			return(-1);
		}

		(*_nup) = (*_nup) + nupi;

		TmpSet1.clear();
		std::set_intersection (OuterBoundary.begin(),OuterBoundary.end(), InnerSubBoundary.begin(),InnerSubBoundary.end(), 
						std::insert_iterator< std::set<int> >( TmpSet1, TmpSet1.begin() ) );
		std::set_intersection (OuterSubBoundary.begin(),OuterSubBoundary.end(), InnerBoundary.begin(),InnerBoundary.end(), 
						std::insert_iterator< std::set<int> >( TmpSet1, TmpSet1.end() ) );

		// Merge the two outer boundaries and remove TmpSet1 from the union

		OuterBoundary.insert(OuterSubBoundary.begin(),OuterSubBoundary.end());	

		TmpSet2.clear();
		std::set_difference (OuterBoundary.begin(),OuterBoundary.end(), TmpSet1.begin(), TmpSet1.end(),
						std::insert_iterator< std::set<int> >( TmpSet2, TmpSet2.begin() ) );

		OuterBoundary.clear();
		OuterBoundary.insert(TmpSet2.begin(),TmpSet2.end());

		// Do the same for inner boundaries

		InnerBoundary.insert(InnerSubBoundary.begin(),InnerSubBoundary.end());

		TmpSet2.clear();
		std::set_difference (InnerBoundary.begin(),InnerBoundary.end(), TmpSet1.begin(), TmpSet1.end(),
						std::insert_iterator< std::set<int> >( TmpSet2, TmpSet2.begin() ) );

		InnerBoundary.clear();
		InnerBoundary.insert(TmpSet2.begin(),TmpSet2.end());

		// Some points in TmpSet1 may still belong to the new inner boundaries
		// 	test if they have a neighbor in the final OuterBoundary

		TmpSet2.clear();

		for (isi=TmpSet1.begin();isi!=TmpSet1.end();isi++)
		{
			(*_irec) = (*isi);
    			(void) GetNeighborhood2 ( &(F->Header) , _irec , TmpSet2 );

			for(isi2=TmpSet2.begin();isi2!=TmpSet2.end();isi2++)
			{
				isi3 = OuterBoundary.find((*isi2));

				if (isi3!=OuterBoundary.end() )
					InnerBoundary.insert((*_irec));

			}
		}
	}

	// Test if we have to increase the value of rigde magnitude for cells in InnerBoundary

	(*_a) = (float)(*_nup) ;
	(void) rec2RC ( &(F->Header) , _root , _nr , _nc );
	(*_d0) = (D->Data)[(*_nr)][(*_nc)];

	for(isi=InnerBoundary.begin();isi!=InnerBoundary.end();isi++)
	{
		(*_irec) = (*isi);
		(void) rec2RC ( &(F->Header) , _irec , _nr , _nc );

		(*_d) = (D->Data)[(*_nr)][(*_nc)] - (*_d0) ;

		if ( (*_d) > (*_k) * sqrt((*_a)) )
			(M->Data)[(*_nr)][(*_nc)] = (*_a) ;
	}

	TmpSet1.clear();
	TmpSet2.clear();
	free(UPSTREAM);
	InnerSubBoundary.clear();
	OuterSubBoundary.clear();

    	return(0);
}


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


int RidgeMagnitudeArea ( int HandleFlowDir , int HandleDistance , int outlet , float k , int HandleOutput ) 
{
	int root , nup, irec, nr , nc ;
	float d0, d, a ;
	FloatRaster * D ;
	FloatRaster * M ;
	FlowDirGrid * F ;
	set<int> InnerBoundary , OuterBoundary, TmpSet1 , TmpSet2 ;
	set<int>::iterator isi, isi2, isi3 ;
	

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

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

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

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


	root = outlet ;

	if ( // RidgeMagnitudeArea_ ( F , D , & k , M , & root , InnerBoundary , OuterBoundary , & nup , TmpSet1 , TmpSet2 )
	     RidgeMagnitudeArea_ ( F , D , & k , M , & root , InnerBoundary , OuterBoundary , & nup , 
				TmpSet1 , TmpSet2 , isi , isi2 , isi3 , & irec , & nr , & nc , & d0 , & d , & a ) )
	{
		SetErreur("RidgeMagnitude : recursive call failed at root cell (%d)\n",root);
		return(-1);
	} 
   
	return(0);
}

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

int RidgeMagnitudeDrop_ ( FlowDirGrid * F , FloatRaster * Z , FloatRaster * M ,
			    int * _root , set<int> & InnerBoundary , set<int> & OuterBoundary , int * _nup ,
			    set<int> & TmpSet1 , set<int> & TmpSet2 ,
			    set<int>::iterator & isi , set<int>::iterator &  isi2, set<int>::iterator &  isi3 ,
			    int * _irec , int * _nr , int * _nc , float * _z0 , float * _z )			
{
	set<int> InnerSubBoundary , OuterSubBoundary ; //, TmpSet1 , TmpSet2 ;
	int * UPSTREAM ;
	int i,nup1,nupi ;//irec,nr,nc;
//	set<int>::iterator isi, isi2, isi3;

	OuterBoundary.clear();
	InnerBoundary.clear();

	InnerBoundary.insert(*_root);
	(*_nup) = 1;

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

	// Extraction of upstream subcatchments
	
	for(i=0;i<nup1;i++)
    	{
		if ( RidgeMagnitudeDrop_ ( F , Z , M , UPSTREAM+i , InnerSubBoundary , OuterSubBoundary , & nupi ,
					TmpSet1 , TmpSet2 , isi , isi2 , isi3 , _irec , _nr , _nc , _z0 , _z ) )
		{
			SetErreur("RidgeMagnitude_ : recursive call failed\n");
			return(-1);
		}

		(*_nup) = (*_nup) + nupi;

		TmpSet1.clear();
		std::set_intersection (OuterBoundary.begin(),OuterBoundary.end(), InnerSubBoundary.begin(),InnerSubBoundary.end(), 
						std::insert_iterator< std::set<int> >( TmpSet1, TmpSet1.begin() ) );
		std::set_intersection (OuterSubBoundary.begin(),OuterSubBoundary.end(), InnerBoundary.begin(),InnerBoundary.end(), 
						std::insert_iterator< std::set<int> >( TmpSet1, TmpSet1.end() ) );

		// Merge the two outer boundaries and remove TmpSet1 from the union

		OuterBoundary.insert(OuterSubBoundary.begin(),OuterSubBoundary.end());	

		TmpSet2.clear();
		std::set_difference (OuterBoundary.begin(),OuterBoundary.end(), TmpSet1.begin(), TmpSet1.end(),
						std::insert_iterator< std::set<int> >( TmpSet2, TmpSet2.begin() ) );

		OuterBoundary.clear();
		OuterBoundary.insert(TmpSet2.begin(),TmpSet2.end());

		// Do the same for inner boundaries

		InnerBoundary.insert(InnerSubBoundary.begin(),InnerSubBoundary.end());

		TmpSet2.clear();
		std::set_difference (InnerBoundary.begin(),InnerBoundary.end(), TmpSet1.begin(), TmpSet1.end(),
						std::insert_iterator< std::set<int> >( TmpSet2, TmpSet2.begin() ) );

		InnerBoundary.clear();
		InnerBoundary.insert(TmpSet2.begin(),TmpSet2.end());

		// Some points in TmpSet1 may still belong to the new inner boundaries
		// 	test if they have a neighbor in the final OuterBoundary

		TmpSet2.clear();

		for (isi=TmpSet1.begin();isi!=TmpSet1.end();isi++)
		{
			(*_irec) = (*isi);
    			(void) GetNeighborhood2 ( &(F->Header) , _irec , TmpSet2 );



			for(isi2=TmpSet2.begin();isi2!=TmpSet2.end();isi2++)
			{
				isi3 = OuterBoundary.find((*isi2));

				if (isi3!=OuterBoundary.end() )
					InnerBoundary.insert((*_irec));

			}
		}

	}

	// Test if we have to increase the value of rigde magnitude for cells in InnerBoundary	
	
	(void) rec2RC ( &(Z->Header) , _root , _nr , _nc );
	(*_z0) = (Z->Data)[(*_nr)][(*_nc)];


	for(isi=InnerBoundary.begin();isi!=InnerBoundary.end();isi++)
	{
		(*_irec) = (*isi);
		(void) rec2RC ( &(Z->Header) , _irec , _nr , _nc );
		(M->Data)[(*_nr)][(*_nc)] = (Z->Data)[(*_nr)][(*_nc)] - (*_z0) ;
	}

	TmpSet1.clear();
	TmpSet2.clear();
	free(UPSTREAM);
	InnerSubBoundary.clear();
	OuterSubBoundary.clear();

    	return(0);
}


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


int RidgeMagnitudeDrop ( int HandleFlowDir , int HandleElevation , int outlet , int HandleOutput ) 
{
	int root , nup, irec, nr , nc ;
	float z0, z ;
	FloatRaster * Z ;
	FloatRaster * M ;
	FlowDirGrid * F ;
	set<int> InnerBoundary , OuterBoundary, TmpSet1 , TmpSet2 ;
	set<int>::iterator isi, isi2, isi3 ;
	

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

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

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

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


	root = outlet ;

	if ( RidgeMagnitudeDrop_ ( F , Z , M , & root , InnerBoundary , OuterBoundary , & nup , 
				TmpSet1 , TmpSet2 , isi , isi2 , isi3 , & irec , & nr , & nc , & z0 , & z ) )
	{
		SetErreur("RidgeMagnitudeDrop : recursive call failed at root cell (%d)\n",root);
		return(-1);
	} 
   
	return(0);
}


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

int LinkBodies_ ( FlowDirGrid * F , FloatRaster * B , FloatRaster * Out , int * _ndown0 , int * _root , char * _existdown , char * _existup , int * _ndown , set<int> & Nodes )
{
	int i,nup1,nr,nc,ndown_i ;
	int * UPSTREAM ;
	char existup_i, existdown_i,b ;

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

	(void) rec2RC ( &(B->Header) , _root , & nr , & nc );
	
	if ( (B->Data[nr][nc]>0 ) | (*_existdown)>0 )	
		existdown_i = 1;
	else
		existdown_i = 0;

	
	// existdown : (down,cell[
	// existup : [cell,up)

	if ( B->Data[nr][nc]>0 )
	{	
		(*_existup) = 1;
		(Out->Data)[nr][nc] = 1;
		b = 1;
		(*_ndown) = 0;
	}	
	else
	{
		(*_existup) = 0;
		(Out->Data)[nr][nc] = 0;
		b = 0;
		(*_ndown) = (*_ndown0) + 1 ;			
	}

	for(i=0;i<nup1;i++)
    	{
		if ( LinkBodies_ ( F , B , Out , _ndown0 , UPSTREAM + i , & existdown_i , & existup_i , & ndown_i , Nodes ) )
		{
			SetErreur("LinkBodies_ : recursive call failed\n");
			return(-1);
		}

		if (existup_i>0)
			(*_existup) = 1;

//		if (i==0 & b==0)
//			(*_ndown) = ndown_i+1 ;		
//		else
//		{
			if(b==0 & (ndown_i+1)<(*_ndown))
				(*_ndown) = ndown_i+1 ; 
//		}
	}

	if ( (*_existup)>0 & (*_existdown)>0 )
		(Out->Data)[nr][nc] = 1;


	if ( (*_ndown)==(*_ndown0) & (*_existdown)==0 )
		Nodes.insert((*_root));

	free(UPSTREAM);

	return(0);
}

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

int LinkBodies ( int HandleFlowDir , int HandleIndic , int outlet , int ndown , int HandleOutput , int ** _NODES , int * _nnodes )
{
	char existup,existdown ;
	FlowDirGrid * F ;
	FloatRaster * B ;
	FloatRaster * Out ;
	int root , ndown0 , ndown_i, nnodes,i ;
	set<int> Nodes ;
	set<int>::iterator isi;


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

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

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

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


	root = outlet ;

	existdown = 0;
	ndown0 = ndown ;
	Nodes.clear();

	if ( LinkBodies_ ( F , B , Out , & ndown0 , & root , & existdown , & existup , & ndown_i , Nodes ) )
	{
		SetErreur("LinkBodies : recursive call failed at root cell\n");
		return(-1);
	}

	nnodes = Nodes.size();
	(*_nnodes) = nnodes ;

	if ( !( (*_NODES) = (int *) malloc( nnodes * sizeof(int) ) ) )
	{
		SetErreur("LinkBodies : could not allocation integer pointer for output\n");
		return(-1);
	}	

	i=0;

	for(isi=Nodes.begin();isi!=Nodes.end();isi++)
	{
		(*_NODES)[i] = (*isi);
		i++;
	}

	Nodes.clear();

	return(0);
}


// Trouver tous les points qui n'ont pas de glace à l'aval et pas de glace à l'amont à une distance de n pixels.




